## This script generates the depictions of the threat and support networks (Appendix Figure 4) for "The Coevolution of Networks of Interstate Support, Interstate Threat and Civil War" 
## This script was adapted from the script in Kinne & Bunte, "Guns or Money," BJPS 2018

library(igraph)
library(sna)
library(reshape2)
library(RColorBrewer)
library(plyr)
library(ggplot2)
library(GGally)
library(grid)
library(gridExtra)

rm(list=ls(all=TRUE))

#dir <- "[Directory Information Here]"
dir <- "C:/Users/rapiduser/Box/Relation creation/JOPreplication/"

setwd(dir)

## Import network/dyadic data
nets.total <- read.csv("relation_network.csv", header=TRUE, row.names=1)

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(nets.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))

## Use last 5 years of the period under analysis
yr <- seq(2006, 2010)

## First sum all support ties during this period
for (jj in yr) {
    net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
    net <- acast(net, ccode1~ccode2, value.var="support_bin")
    out <- emat
    out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
    out[is.na(out)] <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    if ( jj == yr[1] ) {
    	s.net <- out
    } else {
    	s.net <- s.net + out
    }
    rm(net,out)
}; rm(jj)
## Dichotomize to indicate any support during this period
s.net[s.net > 1] <- 1
diag(s.net) <- 0

## Now do the same for threats
for (jj in yr) {
    net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
    net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
    out <- emat
    out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
    out[is.na(out)] <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    if ( jj == yr[1] ) {
    	t.net <- out
    } else {
    	t.net <- t.net + out
    }
    rm(net,out)
}; rm(jj)
t.net[t.net > 1] <- 1
diag(t.net) <- 0

## Measure (in/out)degree in each network
degree <- sna::degree
s.odeg <- degree(s.net, gmode="digraph", cmode="outdegree")
s.ideg <- degree(s.net, gmode="digraph", cmode="indegree")
t.odeg <- degree(t.net, gmode="digraph", cmode="outdegree")
t.ideg <- degree(t.net, gmode="digraph", cmode="indegree")

## Make an edgelist version of each network
s.edg <- melt(s.net); names(s.edg) <- c("ccode1","ccode2","support")
s.edg <- s.edg[s.edg$ccode1 != s.edg$ccode2,] # Drop self dyads
s.edg <- s.edg[s.edg$support != 0, ]

t.edg <- melt(t.net); names(t.edg) <- c("ccode1","ccode2","threat")
t.edg <- t.edg[t.edg$ccode1 != t.edg$ccode2,] # Drop self dyads
t.edg <- t.edg[t.edg$threat != 0, ]

## Code regions based on COW codes, for three axes
reg <- vector(length=length(id))
for (z in 1:length(reg)) {
    if ( as.numeric(colnames(s.net)[z]) < 400 ) reg[z] <- 1
    if ( as.numeric(colnames(s.net)[z]) > 399 &
         as.numeric(colnames(s.net)[z]) < 700 ) reg[z] <- 2
    if ( as.numeric(colnames(s.net)[z]) > 699 ) reg[z] <- 3
}; rm(z)

## Put all the nodal data in a single data frame
node.list <- data.frame(
    name = rownames(s.net),
    s.odeg = s.odeg,
    s.ideg = s.ideg,
    s.deg = s.odeg + s.ideg,
    t.odeg = t.odeg,
    t.ideg = t.ideg,
    t.deg = t.odeg + t.ideg,
    region=reg
)
rm(s.ideg, s.odeg, t.ideg, t.odeg, reg)

## Now calculate two-paths in each network
s.2paths <- s.net %*% s.net; isSymmetric(s.2paths)
t.2paths <- t.net %*% t.net; isSymmetric(t.2paths)

## Function for mapping support 2-path data to edgelist
F1.map <- function(x) {
    data.frame(
        s.2paths = s.2paths[which(rownames(s.2paths) ==
                                      as.character(x$ccode1)),
                                which(colnames(s.2paths) ==
                                      as.character(x$ccode2))]
    )
} 
s.edat <- ddply(
    s.edg, .variables=c("ccode1", "ccode2", "support"),
    function(x) data.frame(F1.map(x))
)
rm(s.2paths, s.edg, F1.map)

## Do the same for threat
F2.map <- function(x) {
    data.frame(
        t.2paths = t.2paths[which(rownames(t.2paths) ==
                                        as.character(x$ccode1)),
                                  which(colnames(t.2paths) ==
                                        as.character(x$ccode2))]
    )
} 
t.edat <- ddply(
    t.edg, .variables=c("ccode1", "ccode2", "threat"),
    function(x) data.frame(F2.map(x))
)
rm(t.2paths, t.edg, F2.map)

## Calculate node size as interpolated degree, bound beteen 0.5 and 1.5
approxVals <- approx(c(0.5,1.5), n=length(unique(node.list$s.deg)))
node.list$s.size <- sapply(
    node.list$s.deg, function(x)
        approxVals$y[which(sort(unique((
            node.list$s.odeg + node.list$s.ideg
        ))) == x)]
)
rm(approxVals)

## Do the same for threat
approxVals <- approx(c(0.5,1.5), n=length(unique(node.list$t.deg)))
node.list$t.size <- sapply(
    node.list$t.deg, function(x)
        approxVals$y[which(sort(unique((
            node.list$t.odeg + node.list$t.ideg
        ))) == x)]
)
rm(approxVals)

## Generate colors for nodes based on degree
F1.col <- colorRampPalette(
    adjustcolor(brewer.pal(9,"YlGnBu")[2:9], alpha.f=0.5),
    bias = 1.5, space = "rgb", interpolate = "linear"
)
colCodes <- F1.col(length(0:max(node.list$s.deg)))
names(colCodes) <- 0:max(node.list$s.deg)
node.list$s.color <- colCodes[ match(node.list$s.deg, names(colCodes)) ]
## Save full spectrum of color codes for legend
s.cols <- colCodes

## Set threat node color based on degree, but using same scale as for support, which
## allows for direct comparisions between networks.
node.list$t.color <- colCodes[ match(node.list$t.deg, names(colCodes)) ]
## Save full spectrum of color codes for legend
t.cols <- colCodes; rm(colCodes,F1.col)

## Assign visual attributes to edges (just support for now but could use 2 paths)
## Use the same scale for both networks, to better
## allow direct comparisons. 
F1.edg <- colorRampPalette(brewer.pal(9,"YlGnBu")[3:9])
colCodes <- F1.edg(length(0:max(s.edat$support)))
names(colCodes) <- 0:max(s.edat$support)
s.edat$color <- colCodes[ match(s.edat$support, names(colCodes)) ]
## Save full spectrum of color codes for legend
s.ecols <- colCodes

## Now the same for threat
t.edat$color <- colCodes[ match(t.edat$threat, names(colCodes)) ]
## Save full spectrum of color codes for legend
t.ecols <- colCodes

rm(colCodes,F1.edg)




################################################################
## Next plot traditional network graphs and degree distributions
################################################################

## Make the support matrix a network object
net.use <- as.network(s.net, directed=TRUE,
                      hyper=FALSE, matrix.type="adjacency")

## Use GGally package to plot ggplot versions of networks
gg.s.net <- ggnet2(
    net.use, edge.color="gray85",
    edge.size=0.1, color=node.list$s.color,
    size=sqrt(node.list$s.deg)/2, size.min=0.001, mode="kamadakawai", alpha=1) +
    guides(color = "none", size = "none") +
    theme(plot.title = element_text(size=12, hjust = 0.03, face="bold"))

gg.s.odeg <- ggplot(as.data.frame(node.list$s.odeg),
                       aes(x=(node.list$s.odeg))) +
    geom_histogram(aes(y=..density..),
                   binwidth=0.5, color="gray25", fill="gray35") +
    xlab("Nodal outdegree") + ylab("Density") +
    theme(
        plot.margin = unit(c(0,0,0,0), "cm"),
        axis.title=element_text(size=10)
    )

a1 <- gg.s.net
b1 <- ggplotGrob(gg.s.odeg)
c1 <- a1 + annotation_custom(grob = b1, xmin=0, xmax=0.4, ymin=0,
                             ymax=100) +
    ggtitle("(a) Support network and\ndegree distribution")
rm(a1,b1)

## Do the same for the threat network
rm(net.use)
net.use <- as.network(t.net, directed=TRUE,
                      hyper=FALSE, matrix.type="adjacency")

gg.t.net <- ggnet2(
    net.use, edge.color="gray85",
    edge.size=0.1, color=node.list$t.color,
    size=sqrt(node.list$t.deg)/2, size.min=0.001, mode="kamadakawai", alpha=1) +
    guides(color = "none", size = "none") +
    theme(plot.title = element_text(size=12, hjust = 0.03, face="bold"))

gg.t.odeg <- ggplot(as.data.frame(node.list$t.odeg),
                       aes(x=(node.list$t.odeg))) +
    geom_histogram(aes(y=..density..),
                   binwidth=0.5, color="gray25", fill="gray35") +
    xlab("Nodal outdegree") + ylab("Density") +
    theme(
        plot.margin = unit(c(0,0,0,0), "cm"),
        axis.title=element_text(size=10)
    )

a2 <- gg.t.net
b2 <- ggplotGrob(gg.t.odeg)
c2 <- a2 + annotation_custom(grob = b2, xmin=0, xmax=0.4, ymin=0,
                             ymax=10) +
    ggtitle("(a) Threat network and\ndegree distribution")
rm(a2,b2)



###################################################
## Print figures
###################################################

setwd(paste(dir, "/Figures", sep=""))

pdf("support_net.pdf",width=11,height=8,paper="special")
par(mar=c(0,0,0,0))
plot(gg.s.net)
dev.off()

pdf("threat_net.pdf",width=11,height=8,paper="special")
par(mar=c(0,0,0,0))
plot(gg.t.net)
dev.off()
